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Abstract 

A new numerical framework, based on the use of a simple first order strongly hyper- 
bolic evolution equations, is introduced and tested in case of 4-dimensional spherically 
symmetric gravitating systems. The analytic setup is chosen such that our numerical 
method is capable to follow the time evolution even after the appearance of trapped 
surfaces, more importantly, until the true physical singularities are reached. Using this 
framework, the gravitational collapse of various gravity-matter systems are investigated, 
with distinguished attention to the evolution in trapped regions. It is justified that in 
advance to the formation of these curvature singularities, trapped regions develop in all 
cases, thereby supporting the validity of the weak cosmic censor hypothesis of Penrose. 
Various upper bounds on the rate of blow-up of the Ricci and Kretschmann scalars and 
the Misner-Sharp mass are provided. In spite of the unboundedness of the Ricci scalar, 
the Einstein-Hilbert action was found to remain finite in all the investigated cases. In 
addition, important conceptual issues related to the phenomenon of topology changes 
are also discussed. 
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1 Introduction 



The diffeomorphism invariance of Einstein's theory of gravity is in an intimate relation with 
the fact that there is a significant redundancy in the representation of the true physical de- 
grees of freedom. It can therefore be a great challenge to carry out a faithful investigation of 
dynamical processes even in case of spherically symmetric spacetimes, in spite of the simpli- 
fications offered by the symmetries. Correspondingly, the selection of the most appropriate 
variables, done by applying a suitable gauge fixing, i.e., the most effective framework to 
carry out the study of a given dynamical system, is considered to be a sort of art. 

In the numerical investigations of spherically symmetric dynamical systems, the method 
of Choptuik — that had been applied first by him in [6] while exploring the critical phe- 
nomenon in the gravitational collapse of various gravity-matter systems — turned to be the 
most successful in the sense that it is still widely used. However, Choptuik's choice has 
both advantages and disadvantages. Perhaps the most important advantage is that one 
has to solve only two first order partial differential equations (PDEs) for the basic metric 
variables 0, which, along with the matter field equations, determine the full evolution of 
the associated gravity matter system. On the other hand, the following objections may 
also be raised. First of all, the radial coordinate r is chosen so that the area A of the 
SO (3)-invariant 2-spheres are given as A = 4-7rr 2 . As it is well-known, coordinate systems 
of this type are not suitable to follow evolution in regions where "trapped surfaces" are 
formed because a coordinate singularity also develops (for a short discussion see, e.g., the 
last paragraph of [53]). 

Another objection is that one of the first order PDEs is hyperbolic while the other is 
elliptic. This means that one of the metric equations — which is a constraint equation — has 
to be integrated on succeeding time level surfaces repeatedly. This process slows down 
time integration and makes it hard to apply the powerful tool of adaptive mesh refinement 
(AMR) which ensures high numerical accuracy in strongly dynamical processes J3j [T2"l [13] . 
There are two side-remarks in order. First of all, although it is possible to implement a 
variant of AMR for the numerical integration of mixed hyperbolic and elliptic equations [33] , 
the precision is partly lost because the applicability of AMR necessitates the extrapolation 
of some variables in time, something which is better to be avoided in a time integration 
process. Thereby, the use of a fully hyperbolic system is preferable in case of numerical 
integration of the field equations based on a finite difference schema. Secondly, it was shown 
by one of the present authors in [33] that by making use of the Kodama vector field as the 
time evolutions vector field, the mixed elliptic-hyperbolic system can be replaced by a fully 
hyperbolic one. However, we have found the inevitable formation of a coordinate singularity 
in numerical simulations, the appearance of which was anticipated in the discussion below 
Equation (4.15) of [33]. The above findings motivate the search for a better analytic set-up 

In this particular case, the basic metric variables are the smooth functions A and B in terms of which 
the spacetime metric can be given as 
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with a more appropriate choice of the basic variables. 

Before proceeding and presenting our proposal for such a choice we would like to mention 
that several attempts have also been made to use reduced versions of the ADM and BSSN 
[TJ 00] formalisms in spherical symmetry, see [5j 122"! [23] . These models have been yielded by 
the reduction of a more complicated evolutionary system, therefore they do not optimally 
fit the spherically symmetric setup. Either they apply coordinates that are not suitable 
to cover regions with trapped surfaces, as it happens, e.g., in [23(23], or they are simply 
too complicated — see, e.g., the basic set of field equations (9a)-(9f) and (10a)-(10c) in [5]. 
They may also suffer from numerical instabilities at the origin because negative powers of 
r appear in the evolutionary equations. 

Thereby, it is of considerable interest to single out a simple and general enough frame- 
work within which time evolution can be investigated on the largest possible part of the 
physical spacetime up to the appearance of true geometrical singularities. To match this 
requirement, we shall start by choosing a fully hyperbolic evolutionary system that is auto- 
matically applicable to describe the evolution in trapped region(s). This choice should be 
such that the equations are free from the numerical instabilities that used to appear in the 
origin in spherically symmetric spacetimes. 

An analytic framework fitting the above outlined expectations may be chosen as follows. 
Based on the results of earlier investigations in spherically symmetric (see, e.g., |29[ 19] [15]) 
and also in generic (see, e.g., [31]) dynamical spacetimes, the metric of the four dimensional 
spherically symmetric spacetime (M,g a b) will be assumed to possess the form 

ds 2 = a f3 2 dr 2 -a dp 2 -r 2 (M 2 + sm 2 $ dip 2 ) , (2) 

where the coordinates r and p label the points of the two-dimensional timelike surfaces 
transverse to the transitivity surfaces of the rotation group, and a, (3, r are smooth functions 
of (r, p). Note that p generally differs from the area-radial coordinate. This condition, as 
we have already mentioned, is necessary to extend the domain of time evolution to include 
trapped surfaces — when they exist. 

By making use of this geometrical framework, gravitational collapse have been investi- 
gated in some simple gravity-matter systems. In the simplest possible case of asymptotically 
flat configurations, the associated time evolution of the system is qualitatively expected 
[51 EJ [TU1 H21 IS] to be as indicated on Fig.Q] where the event horizorH, and the apparent 
horizon are also shown. The latter is foliated by marginally trapped surfaces and it is rep- 
resented by a curve connecting the two ends of the zigzag line depicting the singularity. A 
portion of matter either falls into the singularity or reaches future timelike or null infinity, 
i + or ,y + . The asymptotic structure of the spacetime is expected to approach that of the 
Schwarzschild solution as we get closer to i + along the null generators of [9]. If the 
collapsing matter has no considerable radiative degrees of freedom, the mass of the devel- 
oping black hole is expected to be close to the total mass on the initial data surface. Note 
that spherical symmetry makes it possible to define mass and energy inside an invariant 

2 The event horizon is the boundary of the causal past of future null infinity, Jtf? = <9J _ [J^ + ]. 
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Figure 1: A typical spacetime diagram representing the gravitational collapse of a simple spherically 
symmetric gravity-matter system. In our numerical investigations, we shall focus on the characterisation of 
a gravitational collapse by monitoring the intersections of the Cauchy surfaces and the apparent horizon — 
indicated by circles. 

metric sphere (see Section^]). If the mass inside the marginally outer trapped surface is 
found to approach the total mass while moving outwards along the apparent horizon, then 
we have a strong indication that our spacetime grid covers the truly dynamical part of the 
collapse. Moreover, as it will be demonstrated in Section[71 the lapse function j3 may always 
be chosen such that the Cauchy surfaces can get arbitrarily close to the singularity. 

We were also interested in investigating the time evolution of more exotic initial data 
specifications. Likewise in the standard Friedman- Robertson- Walker cosmological models — 
these are known to be spherically symmetric around any of their spacetime events — there 
is a freedom in choosing the topology of the initial data surfaces. 

The base manifold M of the investigated spacetimes coincides with the future Cauchy 
development of some three-dimensional achronal hypersurface X. Thereby, M = 
and it possesses the product space structure £ x R + . Since the spacetime is spherically 
symmetric, the r = const time level surfaces — these are diffeomorphic to £ — can also 
be foliated by the transitivity surfaces of the rotation group. In virtue of the particular 
form of the applied line element ([2]) the metric induced on the transitivity surfaces of 
the rotation group is given as ds 2 \^> = r 2 (di9 2 + sin 2 ?? dip 2 ) . Thereby, the vanishing of 
r, which is, in fact, the area-radius function, is directly related to the existence of an 
origin. To see how many origins we might have, let us consider some simple choices for the 
topology of the initial data surface S. Whenever E is a connected geodesically complete 
spacelike hypersurface possessing a trivial product bundle structure, its topology is either 
M 3 = [0,oo) x § 2 , S 3 = [0,tt] x§ 2 ,lx § 2 or S 1 x S 2 (see as an illustration the bottom 
line of Fig. [2]). Accordingly, there may be one origin, two origins or no origin at all on our 
initial hypersurface. During time evolution, the geometrical properties of the time level 
surfaces may change such that new origins are produced, possibly indicating a change in 
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Figure 2: Topology changes of the time level surfaces are illustrated by the vertical figure sequences, with 
the suppression of one space dimension. Time progresses from bottom to top. Transitivity surfaces of the 
rotation group are represented by horizontal circles in the first two sequences (starting from R 3 and S 3 ) and 
by vertical circles in the last two sequences (R x § 2 and S 1 x S 2 ). 

the topology, as illustrated on Fig.[2j 

At this point it is important to clear up the some of the related notions and potential 
misconceptions. While in the above sentences the terminology of topology change have 
been used, it should be kept in mind that it is the topology of the limit of time level 
surfaces — comprising the future boundary of the pertinent "domain of dependence" — what 
may differ from that of the time level surfaces. A careful approach is needed in applying 
notions such as domain of dependence, Cauchy development and Cauchy surface which were 
introduced in [25] by Geroch. We would like to emphasise first that there is a significant 
distinction between the notions "initial data surface" and "Cauchy surface" . We may chose 
an arbitrary achronal hypersurface E as an initial data surface. The initial data together 
with the field equations can be used to determine the solution everywhere in the domain 
of dependence D[E]. As opposed to this, the use of the notion of "Cauchy surface" tacitly 
requires additional knowledge about the global properties of the underlying spacetime. In 
particular, see Theorem 11 of [25], a spacetime is known to be globally hyperbolic if and 
only if it possesses a Cauchy surface. 

The confusion might arise — and, in fact, apparently it does arise in many circumstances — 
in consequence of the fact that the domain of dependence D[E] of any initial data surface 
E, is a globally hyperbolic spacetime itself on its own right, i.e., E is a Cauchy surface for 
D[S]. In many cases it is obvious to recognise that there exists a globally hyperbolic global 
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extension^] (M', g' ab ) of (D[T],g a t,) such that the boundary of D[T] is not empty in M', while 
the Cauchy horizon, H[T,] = H + [T,] U -ff~[E] of D[Y,] is obviously empty with respect to 
!?[£]. However, a generic method that could be applied in all the possible cases does not 
exist. Although the concept of "maximal Cauchy development" was introduced in [7], it is 
based on Zorn's lemma, which makes its use difficult in practice. 

As a simple example, consider the maximal analytic extension of the Schwarzschild 
spacetime. In Kruskal coordinates (T, R, ip) the line element is given by 

ds 2 = a dT 2 - a dR 2 - r 2 (dtf 2 + sin 2 tf dip 2 ) , 

32M 3 e~ r / 2GM 

a = , 

r 

where the Schwarzschild coordinates r and t are determined by the following implicit rela- 
tions (see, e.g., [S]): 

(l-—) e r / 2GM = T 2 -R 2 , 
V 2MJ 

-4- = 2 tanh" 1 ( ) 
2M \R J 

By introducing suitable new coordinates r = t(T,R) and p = p(T,R), the line element of 
the Schwarzschild spacetime preserves the form ([2]). The time slicings — two of which are 
indicated on Fig. [3J — and the corresponding lapse functions vary accordingly. In the first 
case (left panel), two points on the S T * hypersurface hit the r = Schwarzschild singularity. 
The appearance of these points can be interpreted as the formation of two origins, since 
r > in all other points of S r * . Correspondingly, the topology E x S 2 of S T0 is replaced by 
the disjoint union of M 3 = M + x S 2 , § 3 and M. 3 at r* (see the third column of Fig. [2]). It is 
not hard to modify S T0 by inserting suitably located further "steps" such that the number 
of §> 3 components of the pertinent £ r * may take an arbitrary integer value. 

Although the time level surfaces on the left panel do not cover the future domain of 
dependence, D + [£ T0 ], of E T0 , it is worthwhile to keep in mind that S T0 is a Cauchy surface 
for the Schwarzschild spacetime. As it is indicated on the right panel of Fig. [31 the entire 
of -D + [£ ro ] can be covered by time level surfaces if the new coordinates r and p, and in 
turn, the lapse /? are chosen properly. There is no topology change and the section of the 
Cauchy surfaces within the black hole region uniformly converge to the r = singularity. 
In section[71 the application of this type of dynamically determined (3 will be demonstrated 
via various examples. 

We would like to mention that the type of topology changes indicated by the left panel of 
Fig. [3J showed up in many of our numerical simulations. We have found that, in consequence 
of the developing inhomogeneities, certain parts of the time level surfaces may get closer to 
the singularity much faster than others. Nevertheless, as shown by the above example, the 
change of the topology refers to the limit of the applied time level surfaces rather than the 
underlying spacetimes. As it will be demonstrated in the succeeding sections, the j3 function 

3 For the precise notion of 'local' and 'global' extensions of spacetimes see, e.g., |36l 137] . 
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Figure 3: Time slicings with and without topology change. On the left, new origins appear (full circles) 
during the evolution. The topology of E T * differs from that of the E T hyp ersurf aces. These time level 
surfaces cover only a part of the future domain of dependence D + [E T0 ]. On the right, it is indicated that the 
entire of D + [E T0 ] can be covered by a properly chosen time slicing such that no topology change happens. 

can be chosen such that the time evolution is slowed down in spacetime regions close to 
the singularity. This way we were able to enlarge the domain of dependences significantly 
and to get arbitrarily close to the developing spacetime singularities in our simulations. As 
a consequence, we were able to investigate the rate of blowing up of the curvature while 
approaching the singularity. 

The structure of this paper is the following. In Section^ with the help of a simple- 
minded approach, we demonstrate the main technical issues in solving the Einstein-scalar 
field equations numerically. These difficulties are resolved in Section[3l where the Misner- 
Sharp mass is introduced as an auxiliary variable to stabilise the applied numerical represen- 
tation near the origin. Here we also describe the "trick" which makes our framework capable 
to follow the time evolution within the trapped region. Section^] outlines the analytic pro- 
cedure that guarantees the required stability of the time evolution in the neighbourhood of 
an origin. In Section[5l a method is presented which provides the opportunity to monitor 
the energy transfer processes. Trapped surfaces, trapped regions and apparent horizons 
are defined in Section[6j Dynamical investigations of gravitational collapse are presented in 
simple and in certain less obvious circumstances are presented in Section[71 The results of 
the convergence and accuracy tests are also presented in this section. The paper is closed by 
summarising our results and by providing a short discussions concerning the interpretation 
and some of the immediate consequences of our findings. 

2 A simple-minded approach 

To have a truly dynamical spacetime, some matter fields also have to be specified in addition 
to the already introduced geometrical framework. For simplicity, we shall restrict our 
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considerations to a real self-interacting scalar field with Lagrangian 

2 = \g ef Ve^f^-V, (7) 

where the potential possesses the form relevant to the Klein-Gordon field, V{ip) = \[P"i\P' 
with mass parameter jjl > 0. Note, however, that V could be chosen to be any sufficiently 
regular function of the tjj field. 

The field equations relevant for the gravity-matter system can be given as 

E a b = Rab ~ ^9abR - 8nT ab = , (8) 

V e V^ + ^ = 0, (9) 



where the energy-momentum tensor reads as 



T ab = VaV'VfeV - 9a 



(10) 



It is straightforward to verify, by a direct calculation, that our basic variables a, (5, r and 
ip satisfy the evolution equations 

a 2 - 8 2 a 2 

o T a T = (3 Opa. p — 87ra(ip T — (3 ip p ) + 



a 

a(3 2 +r 2 - P 2 r 2 (3 T a T + (3 2 (3 P 



+ 2a - ' - ,p + ^" r ^ ^ ap +2af3dp(3 p , (11) 

9 9 U0 2 + r 2 - f3 2 r 2 a 3 T r T , N 

d T r T = (3 2 d p r p + Mraf3 2 V -— — — ^ + —f- , 12 

r p 



d T ip T = f3 A d p ip p -2 — H ~ &tp ' ^ ' 



where the abbreviating notation f a = d a f is applied (for f = a, f3, r and ip). Similarly, the 
constraint equations are 

E T p = 0, E\ = 0, (14) 

where 



pT _ 2 P (d T r p ) -r T a p - a T r p 2r T f3 p if; p ijj T 
h »- ^Wr 



ET = 2a(3 2 (dprp)-r T a T -(3 2 rpap r 2 - r p 2 (3 2 + a(3 2 | /[ ^ 2 + (5 2 ^ p 2 + 2ap 2 V 
a 2 (3 2 r a(3 2 r 2 a(5 2 

(16) 

The number of independent equations is in accordance with the fact that the evolution 
equation (|13j) for ip is simply the actual form of ([9]), and that in a generic four-dimensional 
spherically symmetric spacetime there may only be four independent equations derived 
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from the Einstein's equation ([5)). In particular, equations (|lip and (|12p can be seen to be 
equivalent to the combinations E^ + i? 1 ^ — i? 7 ^ — E p p = and E T T + = 0, whereas the 
constraints are E T p = and E T T = 0. 

Notice that we have no evolution equation for (3. According to the terminology intro- 
duced by Priedrich in [21], this function can be seen to play a role analogous to the "gauge 
source function". It is freely specifiable, provided that it is sufficiently regular on the do- 
main of time evolution. By making use of this freedom, we were able to slow down evolution 
in our numerical simulations at the parts of the time level surfaces whenever they have got 
close to the singularity. 

Defining the coordinate speed of light, cupi = dp/dr along the radial null geodesies in 
the t — p part — requiring the vanishing of the pertinent part of the line element — , it is 
straightforward to see that c\ T>p \ = j3. Thereby, as far as (3 is guaranteed to be bounded so 
is c\ T p], In our simulations, (3 was guaranteed to be less than or equal to 1.2 everywhere. 
It is also remarkable that the system of evolution equations may be viewed as a system of 
coupled non-linear wave equations for the variables a, r and ip. These type of equations are 
known to possess a well-posed initial value problem (see, e.g., [llj). In addition, it can be 
checked by a direct calculation that whenever the evolution equations are satisfied — and, 
in turn, their derivatives also vanish — the following set of first order strongly hyperbolic 
evolution equations can be derived 

2 r , (otP + ad T )E T n - (2a /3 + a(3 p )E T T 
9tE p - d p E\ + - [r T ET p - r p E\\ + ^ Pt) p J pP ^L^L = 0, (17) 

o 2 r on a T E T r — (3{a (3 + 3af3 )E T n 

d T E\ - d 2 d p E\ + - [r T E\ - f3 2 r p E\] + ^— — = . 18 

t H a 

Since these equations are linear and homogeneous in the variables E T T and E T p , they ensure 
that the constraints propagate with the evolution, therefore it is enough to impose them 
only on the initial data hypersurface. 

In spite of the attractive features of the above outlined evolutionary system — see equa- 
tions (|imi3p and (fT4|) along with (fT7|) - ([T8|) — , it has been found to be unstable in numerical 
simulation because truncation errors grow rapidly near the origin(s), i.e., where r tends to 
zero. The instability is caused by terms of the form 

- — ^r^> (19) 

appearing in (|12p with n = 1 and in (jlip with n = 2. The nominator of (|19p consists of 
three independent variables, each of which contributes to its numerical error. Although it 
should vanish rapidly while approaching the origin, it does not exactly do that in practice 
and, in turn, its slightest error is amplified enormously via the division by r n . 

Increasing the resolution does not help in such a circumstance, the simulations crash even 
earlier. We tried to cure this instability in various ways such as extrapolation from outer 
points, decreasing the order of finite differences near the origin, increasing the numerical 
dissipation term, and also by using different parametrisations of the metric which were 
supposed to behave better close to the critical points. Neither of these attempts turned out 
to be successful. 
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3 Stabilising with the Misner- Sharp mass 



By the inspection of the term (|19p responsible for the above mentioned instabilities and of 
the Misner-Sharp mass m [30] defined by the relation 



m = r - (l + g ab d a rd b r) = ( ~" ' ■ T - " ■ p ] , (20) 



af3 2 + r 2 - j3 2 rj 
2 I a^ 2 

it seems to be rewarding to introduce m as an auxiliary variable. By deriving it with 
respect to r and p and by using the pertinent form of Einstein's equations (JSJ) to eliminate 
the second derivatives of r, an evolution and a constraint equation can be deduced: 



d T m = -{r T EP p - r p E" T + 8vr (r> T% - r p T" T ) } , (21) 
d p m = j {r p E\ - r T E? p + 8tt (r p T\ - r T T T p )} . (22) 



What makes the introduction of m even more preferable is related to the following obser- 
vations. First of all, Einstein's equations E ab = can be seen to hold whenever 

d T a T = P 2 d p a p + 4m ^ 2 " 2 + 8vra 2 /3 2 (t% + T% - T\ - T% 



a 2 T -P 2 a 2 p [3 T g T + p 2 p p ^ 
H H ^ + 2a pd p p p , (23) 



a f3 

2m (5 2 a f3 T r T + (3 2 P p r p 



d T r T = (3%r p + AvraP 2 (T; + T? p ) - ^1 + ^IlZ^ELP , (24) 

d T m = 4vrr 2 (r T T p p - r p T P T ) , (25) 
d p m = 47Tr 2 (r p T T T -r T T T p ). (26) 

It is important to emphasise that instead of solving the constraint equations E T T = and 
E T p = on the initial data surfaces, it suffices to solve (|26p . provided that the evolution 
equations (p3]) - ([25|) hold. To see this, recall first that (|2"3]) is equivalent to 



E T T + E p p = 0. (27) 

In virtue of {SJ, (J22J), (g5D and H2BJ, we have that 

r T EP p -r p Ef T = 0, r p S T T - r T S T p = 0. (28) 

Equations (|27j) - ()28j) along with the relation _E P T = —f3 2 E T p , imply then that for the variables 
E T T , E T p , the homogeneous linear equations 

v T E\ - (3 2 r p E T p = 0, r p E\-r T E T p = (29) 

hold. The algebraic sub-determinant of this system is r 2 — (3 2 r 2 , which will be shown in 
Section [6] to vanish only at isolated marginally trapped surfaces. Therefore the system ()29p 
possesses only the trivial solution which justifies our claim. 
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Restricting our considerations to the case of a self-interacting scalar field again, the 
components of the energy-momentum tensor can be given as 

K = - 2 \/± + v, (30) 

n = -in = (3i) 



1 V? + f3 2 ^l 

T p = M + V , (32) 



7* = T% = -- YT ^ P +V. (33) 



a(3 2 

1 g ~ 

2 a/3 2 

In this case, the complete set of evolution equations consists of (f23|) - ([25|) and (fl~3"j) . 

Note that the introduction of m makes the use of (|23p superfluous since a can be 
determined from ()20p once the fields m and r are known, However, this way we get a 0/0 
type formula for a at 'marginally trapped surfaces' (for their precise definition see Section[6]) 
and an instability in the numerical simulation, thereby making it impossible to study the 
evolution inside the black hole region. In order to overcome these difficulties, we kept (|23p 
as an evolution equation for a. Then, as it will be demonstrated in Section[7] (see Fig. [7]) we 
could guarantee that our numerical framework does possess the desired capabilities without 
the loss of numerical accuracy anywhere in the domain of dependence. 

By making use of standard procedures (described in details, e.g., in [H]), a first order 
evolutionary system for our spherically symmetric dynamical system can be deduced as 
follows. In addition to m, a, r, and ijj, the first derivatives a T , a p , r T , r p , ip T and tp p are also 
considered as being independent variables, and the evolution equations are complemented 
by the relations 

(34) 
(35) 
(36) 
(37) 
(38) 
(39) 

Note that the last three relations are, in fact, the integrability conditions for a, r, and ip. 
As a byproduct of this reduction process, in addition to the true dynamical constraint ([26]) . 
we also have to take care of the trivial ones 

d p a = a p , (40) 

d p r = r p , (41) 

d p i; = V (42) 
The yielded first order system can then be put into the form 

d T u = Ad p u + B, (43) 





= a T , 


d T r 


= r T , 


d T ijj 


= A, 


d T a p 


= d p a T , 


d T r p 


= d p r T , 


d T ip p 


= d p tp T . 
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where the associated 10-dimensional vector variable u, the 10 x 10 matrix A and the 10- 
dimensional source vector B are given as 



/ m\ 
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(44) 



B 



2nr 2 
"off 1 



«r-ff 2 Qp , a T Pr+P 2 a P P P 



+ 



P 



+ 2af3 8 P P P - 8ira(3 2 (tf - + 




_2jn^ + p T r T +p*f} p r p + y 





^ T r T -P 2 j] p r p (3 T i>T+P 2 P P i} p 



V 



p 



o 



(45) 



Since the eigenvectors of the matrix A comprise a complete system and its eigenvalues 
are all real, this first order system is strongly hyperbolic [26] @ It is known that both the 
analytic and numerical well-posedness of the initial value problem is guaranteed for this 
type of evolution equations [351 ES]- An additional preferable property of the investigated 
dynamical system is that for the coordinate speed of light, cupi < /3 max holds, where /3 max 
denotes the maximum value of f3. 



4 Stability in the origin(s) 

The behaviour of the basic field variables m, a, (5, r, ip close to an origin located, say, at 
p = po can be explored by substituting the form of Taylor expansion 

P) = E h/k(r) (p - Po f + O f ((p - P0 ) m ), (46) 

k=0 

4 Note that defining & p = f3ot p , f p = (3r p and -ijj p = (3ip p , the above set of field equations could 
be put into the form of a first order symmetric hyperbolic system for the vector valued variable 

(m, a, a T ,a p ,r, r T ,f p , ip, ^ T ,ip p ) T . 
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into the field equations for the functions m, a, ft, r, ip, and also by requiring the solutions to 
be at least of class C e for some IgN. Assuming that the solutions are at least of class C 4 
in a neighbourhood of po, the resulting relations are 

m = I m 3 (p - Po f + O m {{p - Po f) (47) 

a = a + - a 2 (p - po) 2 + i^a 4 {p- p ) A + O a ((p - p ) 6 ) (48) 

ft = Po + \fti (P - Po? + y A Pa (P - Po) 4 + O p {{p - Po f) (49) 

r = sgn(p - p )^/a~o~ (p - po) + - r 3 (p - p ) 3 + O r {(p - p ) 5 ) (50) 

^ = ih + \ih{p- Po) 2 + ^4(P- Po) 4 + Oi,((p ~ Pof) ■ (51) 

By making use of the indicated parity properties of the functions a, ft, r and ■0, the grid 
boundary at the origin can be treated exactly the same way as it was done in \19\ I20j . 
In particular, the method described in details in Section D of [20J makes it possible to 
determine all the spatial derivatives in a neighbourhood of an origin, by making use of a 
symmetric stencil. 

Note that the factor sgn(/> — po)\A*o m the first order term in (|50p has a simple geomet- 
rical explanation. In order to have a regular origin, i.e., to avoid the appearance of conical 
singularities, we have to guarantee that while approaching to the origin, p — > po, the ratio of 
the circumference of an infinitesimal, origin-centred circle and of the corresponding proper 
radius tends to 2tt, i.e., 

2nr 2ird n r 
i rp i — i I - ? P —~^ - 2vr, 52 

where l'Hospital's rule is also applied. 

Having the relations (|47[) and (|50p . it is also straightforward to see why the system 
introduced in Section[3] has been found stable in our numerical simulations. One of the 
reasons for this unexpected stability is related to the fact that the term f)19|) — which was 
responsible for the numerical instabilities in case of the simple-minded system — reads as 



aft 2 + r 2 T - ft 2 r 2 p 2 m aft 2 m 3 aft 2 



r n+l 



"*3"^ | |2-n /ro\ 

Sa^+^Z 2 ~ ■ (53) 



Since n takes the values 1 and 2, the term (|19p is replaced by a completely regular one in 
our new setup. It is also important to note that the Misner-Sharp mass m is subject to a 
first order PDE (|25p which does not involve a singular term at all. 

If there is an origin, it is useful to have the limiting form of the source vector B appearing 
in (|43p . By making use of (|47p -(|5ip. it is straightforward to verify that the non-trivial 
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components B aT , B Tt and B^ T of B, as given in (|45p . tend to the following regular limits 



3 (d T a ) 2 (d T a )(d T p ) 
p^po " T 2 «o A) 



hm B aT = - 1 - 2sgn(p - pojV^oA) r 3 + «2 /?o 



+ 2a (3 f3 2 - 8ira (d T ^ ) 2 , (54) 
lim B r = , (55) 

P~*P0 

hm R/^ = h 2/? V>2 - aoPo ~5T • (56) 

5 Energy balance 

As it is well-known, energy transfer processes cannot be investigated in general relativity, 
since appropriate quasi-local quantities have not been found as yet. However, this can be 
done in spherically symmetric dynamical spacetimes, using the Kodama vector field [29] 

K a = e ea d e r, (57) 

where e ab = q ae q b f e e f and e a b denotes the volume form associated with the metric q a b 
induced on the two-dimensional surface transverse to the 50(3) group orbits. It was shown 
in [29] that K a is divergence free, i.e., V e K e = 0, and also that G e *V e Kj vanishes, where 
G a b denotes the Einstein tensor R a t> — ^Rg a b. It follows then that the vector field 

ja = T a bK b ( 5 g) 

is also divergence free. Accordingly, even though our spherically symmetric spacetime is 
fully dynamical, J a behaves as a locally conserved energy flux type vector field. Then, by 
making use of Stokes' theorem and by choosing Q to be a four dimensional spacetime region 
with boundary dQ, we have that 

V a J a e = [ n a J a i = 0, (59) 
Q JdQ 

where e denotes the 4-volume element while e is the 3-volume element induced on the 
boundary dQ of Q, which can be given as i a b c = £ e abc ne , where n a is the (outward pointing) 
unit normal 1-form field on dQ. 

With our particular choice of coordinates, the components of the Kodama vector field 
can be given as 

K a = (\ ~— R , 0, o) , (60) 

which, in virtue of (|58[) . implies that the non-zero components of the divergence free vector 
field J a read as 

J T = ^ [T\r p - T T pT T j and J? = JL [ T P T r p - T" p r T ] . (61) 
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Figure 4: The boundary dQ = dQ T2 U c)Q P2 U dQ T1 U 9Qpi of the shaded spacetime domain, enclosed by 
the t = Ti, r = T2, p = pi and p = p2 hyp ersurf aces. 



To characterise the energy transfer processes in a spacetime domain Q bounded by 
sections of r = const and p = const hypersurfaces, as it is depicted on Fig.[H we introduce 
the following auxiliary notations. 

Denote the origin-centred ball of radius p on the r = f hypersurface by B(f,p), and 
the part of the cylindrical hypersurface p = p between r = r\ and r 2 by C(t 2 ,ti, p). Then 
the boundaries dQ n (i = 1,2), as shown on Fig.Hl are represented by the differences 
AB(ri, p 2 , pi) = B(Ti,p 2 ) \B(ri,pi) of the balls Bfa, p 2 ) and B(ri,pi), while dQ Pi are the 
cylinders C(t 2 ,ti, pi) connecting the boundaries of the balls B(t 2 , p%) and B(t%, pi). 

Then the energy contained in AB(f, p 2 , pi) can be given as 

E(f;p 2 , Pl ) = I nPn^, (62) 

JAB(f lP2 ,pi) 

where is the volume element on the r = f hypersurface and na is its unit normal 
1-form field, i.e. 

(63) 



n W = 



>/l(V c r)(V e r)| 

Similarly, the energy transported through the cylindrical hypersurface C(t 2 ,t\,p) is 

5(^,ti; p) = / n®r~eW, (64) 

JC(V2,T1,P) 

where is the volume element on the p = p hypersurface and is its unit normal 
1-form field, i.e., 



(65) 
p=p 
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Using these notations, the energy balance equation (|59p takes the form 



E(t 2 ;P2,Pi)-E(t 1 ;p2 1 Pi) + S(t2,t 1 ;p 2 )-S(t2,t 1 ;p 1 ) = 0. (66) 

The following relations can be verified by straightforward calculation: 

n« = (v^A 0,0,0), n</>) = (0,v^,0,0), (67) 

e ~£ } 7 = V^r 2 S m9-(dp) a A(d$)pA(dip) 7 , (68) 

Va/?r 2 sin(9- (dr) Q A (dtf)^ A (d<p) 7 . (69) 



(p) 



Thus £7(t; P2, Pi) and 5(t2, Ti; p) can be given as 

E(f;p 2 ,pi) = 4tt / (a(3r 2 J T )\ T=f dp, (70) 

■/ 01 



01 

T2 

2 



S{T2,n;p) = 4vr / (q /3r^ JP)\ p=p dr . (71) 



n 



It is important to mention here that, in virtue of the relations (|26|) . (|61|) and 
the Misner-Sharp mass function m can be given on a r = f time level surface as 

m(f,p) = E{f;p, Pl =0) = ?n(f,0)+47r f {a (3 r 2 J T )\ T= r dp . (72) 

J o 

Accordingly, in virtue of (f68|) - ([69j) and (|72p. the expression e GM = ^/aj5J T could also be 
interpreted as the combined energy density of our composed gravity-matter system. 

6 Trapped surfaces 

To define trapped surfaces, start with a smooth orientable 2-dimensional compact manifold 
5? with no boundary in a 4-dimensional spacetime (M,g a b). Let n a be a smooth future 
directed non-vanishing null vector field on which is normal to % i.e., g a bn a X b \y = 
for any vector field X a tangent to S fi . Consider then the null hypersurface M generated by 
geodesies starting on 5? with tangent n a . By parallel propagating n a along these generators, 
it can be extended onto N . Denote by u the synchronised affine parameterisations of these 
null geodesies. The hypersurface M is smooth in a neighbourhood of and it is smoothly 
foliated by the u- level surfaces. Denote by e q the volume element associated to the metric 
q a b induced on these 2-dimensional surfaces. Then the null expansion 8^ with respect to 
n a is defined as 

£ n e q = 9^ e q , (73) 

where £ n denotes the Lie derivative with respect to the null vector field n a . 

It is straightforward to see that the sign of 8^ remains intact under a positive rescaling 
of n a , i.e., whenever n a is replaced by n la = fn a , where / is a sufficiently smooth positive 
function on S fi . Therefore, as far as the sign of 9^ is concerned, n a and n' a may be 
considered to be equivalent. It is well-known that there always exist two equivalence classes 
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of future directed non- vanishing null vector fields on 5? which are normal to 5? . Let n\ 
and n°L be two such vector fields and denote the corresponding null expansions by 9+ and 
According to the original definition of Penrose [32], a 2-dimensional surface y in a 
4-dimensional spacetime is future trapped or untrapped if both of the future directed null 
geodesic congruences orthogonal to 5? are converging, or one of them is diverging while 
the other is converging at 5? , i.e., 9± < or 9 + and 9- have opposite signs. If one of 
the expansions vanishes identically while the other is everywhere non-positive, then 5? is 
called future marginally trapped surface. Past trapped and past marginally trapped surfaces 
can be defined analogously by reversing the time orientation applied tacitly in the above 
definitions. 

Choose now 5? to be a 2-sphere of radius r, invariant under the 50(3) symmetry of our 
spherically symmetric spacetime. In virtue of the p — r part of the metric ([2]), an obvious 
choice for the non- vanishing future directed null vector fields n+ and n a _ on 5? is 

(74) 

while the 2-volume element e q on the 2-spheres ^ T ,p foliating the J\f± null hypersurfaces is 

e q = r 2 sini? • e, (75) 

where e denotes the Levi-Civita symbol. It is straightforward to check that the pertinent 
null expansions are 

9± = - r • (r> ± Pr p ) . (76) 

Equations (|20p and (j76[) imply that the following relation holds on 5?: 

2m 2 r 2 - r 2 r 2 
1-— = P P R2 = - '.O+O.. (77) 
r a/3 2 4 a/3 2 

Since all the functions a, (3 and r 2 are positive, apart from origin(s) or the spacetime 
singularity, an SO (3) invariant 2-sphere 5? of radius r is untrapped or trapped (future or 
past) if and only if 1 — ^ is positive or negative, respectively. 

We shall call a spacetime region future trapped if each point of it belongs to an 
SO(3) invariant trapped 2-sphere. The past boundary, d~ S?+ = 3 '+ n J - [=5^+], of such a 
future trapped region will be referred as future apparent horizon i^+Jfl The notion of past 
trapped region, =^L, and past apparent horizon, =sz/_, may be introduced analogously. 



a Apart from the case of asymptotically flat spacetimes, the '+' and '— ' signs have nothing to do with 
outward and inward pointing directions. Although these directions cannot be defined for generic spacetimes, 
it is worthwhile to mention that an adequate quasi-local notion of outward and inward direction can be 
introduced — without referring to global properties of the underlying spacetimes — provided that attention is 
restricted to untrapped or marginally trapped surfaces (see, e.g. |35]). 

6 It is important to keep in mind that whenever non-invariant topological 2-spheres are involved, the 
extent of the trapped region might be larger than the one determined by S'0(3)-invariant 2-spheres [2]. 
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7 Dynamical investigations 



The first order strongly hyperbolic system described in Section[3] and given in the form (|43p 
was solved numerically with our finite difference code called GridRipper AMR [TJ] (see 
also [1'2[ I13j). The time integration process uses the method of lines based on a fourth 
order Runge-Kutta scheme |26] — for a detailed description of this method, for the case of 
fixed spatial resolutions, see, e.g., [20] — and as it is indicated by its name, adaptive mesh 
refinement techniques are also incorporated. 



7.1 Initial data 

Before determining the time evolution of our system using equation (|43[) . we have to specify 
suitable initial data on a spacelike hypersurface So — which hereafter will always be chosen 
as the r = hypersurface — so that all the constraint equations are satisfied. Recall that 
the components of the vector variable u in (|4"3"]l are m,a,a p ,a T ,r,r p ,r T ,ip,ip p and Vv and, 
as it follows from the discussion in Section[3l we apparently have only the constraint 

27rr 2 

d p m = —p- [r p (^ + + 2a/3 2 V) - 2r r t/v ^ p ) (78) 



beside the three trivial ones (|40p - (|42|) . However, the price to be paid for using the "super- 
fluous" evolution equation f|23|) is the appearance of an additional constraint (|20p . which 
implies that only two of the four variables, a,a T ,r and r T , may be chosen freely. To sat- 
isfy all constraints on the initial hypersurface, we have applied either of the following two 
selection procedures. 

A. Fixing of r and r T from among the functions a, a T ,r and r T on Eo- Start by specifying 
the sufficiently regular but otherwise arbitrary functions r,r T ,ip and tp T , then choose 
r p and ip p so as to satisfy the trivial constraints (f4Tj) and ([4lZ]) . By substituting the 
relation 



2m\- 1 P 2 r 2 p -r 2 
r ) f3 2 



a= 1-— ' > (79) 



into (|78p . the initial data for m can be determined numerically. Then the initial a 
and a p are determined by (|79p and the trivial constraint (|40p . respectively. Finally, 
a T is the r-derivative of (|79"p . 



a T = 2 



(r m T — m r T )(f3 2 r 2 — r 2 ) + r(r — 2m) (/3 2 r p • d p r T — r T ■ d T r T + r 2 f3 T /f3) 

(3 2 (r - 2m) 2 ' 

(80) 

B. Fixing of a, r, ip and ifi T on So- Choose a p ,r p and ip p so as to satisfy the trivial 
constraints (00]) -(g2]). From ([20]), we get 



(81) 



18 



where the sign depends on the specific physical problem we intend to investigate. By 
substituting the right hand side of (|8T|) into (fTH|) . the yielded equation can be solved 
(numerically) for m on the initial data surface So- Then the initial data for r T is 
determined by the relation (|81|) . Finally, a T can be determined similarly as in case A. 

Procedure A. is simpler but B. provides a convenient control of both the temporal geo- 
metrical setup and the energy distribution of the matter field. The latter is also more 
appropriate to investigations of expanding cosmological models. Other choices may also be 
possible, however, the use of these two cases was found to be completely satisfactory. It is 
worth to be emphasised that either of these initial data specifications is suitable to host all 
the possible Einstein-scalar field systems in the selected setup. 

To summarise, we have free control only of the four dynamical variables r,r T ,ip and ip T 
or a, r, if; and ip T in the two cases, respectively. In addition, the specification of (3 is also 
required. Once these functions are specified, they determine all other functions by either of 
the outlined procedures on So- In our numerical simulations, we always used (3 = 1 on So- 

7.2 Dynamical lapse function 

As it was emphasised earlier, the lapse function (3 is freely specifiable. Also, as it was 
indicated by the simple example in the introduction, the extent of the spacetime domain 
covered by the time level surfaces is influenced significantly by the chosen lapse function. 
To achieve the largest possible domain, we use a dynamically determined (3 satisfying the 
evolution equation 

(82) 

where p is a positive real parameter, starting with the initial value /3|s = 1. The value of 
the parameter was chosen to be p = 10 -4 in all of our simulations. 

Equation ([82]) was motivated as follows. To slow down the evolution close to the devel- 
oping singularities, we need a function which tends to zero there rapidly enough. For this 
reason, first we applied (3 = exp(— pm/r 3 ). Besides decaying fast enough, this function is 
also regular at the origin since, in virtue of (|47|) and (|50|) . m/r 3 remains finite while p — > pq. 
The r-derivative of this (3 is "almost" (|82p . but without the factor (r T /r) 3 on the right hand 
side. The inclusion of this factor was found to be necessary in stabilising the time evolution 
close to the singularity. Its use is supported by the observation that in a neighbourhood 
of a point where the function r takes its minimum the factor (r r /r) 3 gets to be sufficiently 
small providing thereby some additional slowing down to the evolution of (3 there. 

7.3 Gravitational collapse with spatial topology M 3 = [0, oo) x § 2 

Let us start by investigating the gravitational collapse of a massive real scalar field in 
the "conventional" case, i.e., with the assumption that the initial data is specified on a 
hypersurface So possessing the topology of M 3 . According to procedure A, we specify two 
metric functions on Eo as 

r{p) = p, r T {p) = 0. (83) 



m T r — 3 m r r 
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The scalar field ip is chosen to be a smooth hunch with compact support 

„ fcexp (d + 7 — h \t ts ) if \p — a\ < b, , , 

[0 if \p-a\ > b, 

with vanishing r-derivative on So 

MP) = 0. (85) 

The self-interaction potential is assumed to possess the form V{ip) = ^p 2 i( 2 . The parameter 
values are 

a = 0, 6 = 70, c = 0.0795, d = 100, fi = 0.8 . (86) 

In the numerical simulation, AMR was used with 7, 000 spatial points on the base grid and 
five refinement levels (refinement ratio 2). 

To demonstrate the significance of the freedom in specifying (3, we shall compare the 
result yielded in two cases, (3 = 1 and the dynamical lapse determined by (|82p . The initial 
state is the same in these two cases, as shown on Fig.0 where the initial matter and gravity- 
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Figure 5: The matter and gravity-matter energy density distributions associated with a shell of radius p 
on the initial data surface Eo. 

matter energy density distributions associated with a shell of radius p, S M = An r 2 yfaT T T 
and t§ GM = 4tt r 2 a f3 J T are plotted. The p integrals of these quantities are, respectively, 
the matter and gravity-matter energy contained in an origin-centred ball of radius p on the 
r = const, time slices, E M (r, 0, p) = J Q P £ M dp and E GM (r, 0, p) = J P & GM dp. 

7.3.1 Time evolution with unit lapse 

The time evolution of this initial state, with (3 = 1, yields the gravitational collapse of the 
scalar field until a spacetime singularity is formed at t* ~ 17.583. The change of the metric 
function r during the pertinent time evolution is depicted on Fig.[6l where the appearance 
of the second origin is clearly visible at p* ~ 2.1. 
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Figure 6: The area radial coordinate function r(p) is plotted on various time slices linearly (left) and 
logarithmically (right). A new origin can be seen to develop at p* w 2.1 on both plots. 

Fig. [71 shows q, a' = (r? — r^)/(l — 2m/r), l — 2m/r and r^ — on an intermediate time 
level surface where trapped surfaces are already present. The value of a is calculated two 
ways: by solving the evolution equation (|23p and by using the algebraic relation (|79p . At 
the first sight, the two calculations yield approximately the same result. However, there is 
a tiny difference close to the marginally trapped surfaces which, in virtue of (|76p and (|77p . 
are associated with the zeros of 1 — 2m /r and r? — in the present case. This difference is 
enlarged by p-derivation, thereby a p and a' p differ significantly here. Whenever we started 
to evolve the system using (|79p . the error became enormous in a couple of time steps, 
making it impossible to carry on with the time evolution in the trapped region. To reach 
long-term stability even in regions with trapped surfaces, it turned out to be essential to 
apply as an evolution equation. 

Fig. [8] is to demonstrate that the combined gravity-matter energy is preserved with high 
precision during the entire history. More precisely, the energy balance relation (|66p . defined 
by making use of the Kodama vector field in the previous section, is found to hold during 
the entire evolution. The wiggly ends of the curves at r ~ r* are associated with the very 
appearance of the scalar curvature singularity at (p*,r*). In the evaluation of the energy 
balance relation (|66p . the particular choices p\ = and p 2 = 69 have been made. Moreover, 
E and Eq are defined as 

E = E(t; p 2 , pi) + S{t, 0; p 2 ), E = E(0; p 2 , pi). (87) 

Notice that doubling the number of grid points shifts the curve downward by about factor 
1/16. Therefore Fig. [8] provides a justification that our numerical method is fourth order 
accurate even with respect to E, which is a complicated function of the elementary field 
variables. Note that the experienced convergence order is in accordance with the applied 
fourth order finite difference schemes, both for the spatial derivatives and for the time 
integration. 

Fig. [9] shows that the constraints ([26]) and ([79]) are preserved with high numerical accu- 
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Figure 7: Functions r 2 p — r%, 1 — 2m/r, a and a' are plotted on a E r hypersurface. On the right, the 
corresponding difference Aa = a — a' is shown. Thin vertical broken lines indicate the location of marginally 
trapped surfaces. Although the error of a' appears to be small, the spikes on the right plot indicate that its 
p-derivative has large errors at the marginally trapped surfaces. 



racy even on the last time slice before reaching the singularity^ The Misner-Sharp mass m, 
its derivative, m p , and the a metric component are also plotted for comparison. Note that 
m appears to blow up at the point where the singularity is about to form, its maximum 
value reaches m max ~ 5 x 10 2 on the time slice shown. The error of the constraints (thick 
solid curves) also have peaks in this point. 



7.3.2 Time evolution with dynamical lapse 

In case of dynamical lapse, the value of the Ricci scalar is of the order 10 10 within the 
interval 3^/9^17atr = 29. Constraint preservation on this time level surface is shown 
by Fig.CEuT 

On Fig.[TT]the time evolution of the gravity-matter energy density distribution associated 
with a shell of radius p, <§ Ghn is shown. This spacetime diagram depicts the variation of 
S GM above the p — r c coordinate plane, where r c denotes the conformal time defined, along 
the constant p-coordinate lines, as r c = J*J~ (3dr. The use of r c is advantageous since, 
as it is indicated by Fig. 1121 the radial null geodesies (see the thin dotted lines) are well 
approximated^! by r c ± p = const lines. On top and bottom, the orthogonal projections of 
the world sheet of $gm can be seen. The curvy boundary represents the singularity. It is 
remarkable that $gm is positive everywhere outside the apparent horizon. However, beyond 
this boundary, it starts to oscillate such that the amplitude of this oscillation is increasing 
towards the intersection of the singularity and the p = line. It is important to note that 
in spite of this oscillatory behaviour of $GMi the Misner-Sharp mass is always nonnegative. 



7 Note that instead of the p-dependent values of the constraint violations, the p-dependence of their 
averaged values are shown, in order to smooth out the noise, on Figs. I9I1T71 The averaging for each p is done 
by using the interval [p — 0.05, p + 0.05] around p. 

8 To see this recall that along radial null geodesies /3dr±dp = 0. This, along with dr c = /3dr+ [J Q T /9 P dr] dp 
and that the r-integral of f3 p is negligible, implies that in the p — r c coordinate plane they satisfy the 
differential equation dr c /dp = =pl + [J^ /3pdr] w =pl. 
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On Fig. 1121 the location of the apparent horizon is plotted on the p-T c coordinate plane. 
In accordance with the generic expectations which are also justified by Figs.[T]and[71 there 
are two marginally trapped surfaces on each r = const time level surface. The p-coordinate 
of the inner one is decreasing while that of the outer one is increasing. The plotted radial 
null geodesies indicate that the apparent horizon on the left panel is everywhere achronal 
as it is expected. It is also remarkable that the outer part of the apparent horizon tends 
to become almost null very rapidly. On the right panel of Fig. [T2] the Misner-Sharp mass is 
plotted as a function of p along the future apparent horizon, As it is clearly visible, m 
is increasing monotonously along and it tends to the value of the ADM mass m ADM = 
m(r = 0, p — > oo) which is equal to the Misner-Sharp mass m(r = 0,p max ) in the present 
case since the matter field is of compact support. As the area of the apparent horizon 
is A = Attt 2 and r/2 = < m ADM , this plot is also in accordance with the Penrose 

inequality \J A/1Qit < m ADM . 

As mentioned above, a true scalar curvature singularity develops by the end of the time 
evolution inside the black hole region when both the Kretschmann scalar, R a b c dR abcd , and 
the Ricci scalar curvature, R 3C = g ab R a b, tend to infinity. As it is justified by the log-log 
plots on left panel of Fig.[13l the blow up rate of R a bcd R abcd , R sc = g ab R a b and m, at fixed 
p = 2.5 location, are r~ 8 , r -4 and r -0,95 , respectively. On the right panel of Fig. IT3l the in- 
dependence of the critical exponents 7 m , and jk — 7 m of the blow up rates m ~ A m r 7m 
and {R a i 1C( iR abcd ) 1 ^' 2 ~ Ak r" /K are shown. It is visible that these curvature scalars, along 
with the Misner-sharp mass, diverge faster close to the origin. We would like to point 
to the fact that the blow-up of the Kretschmann scalar is completely consistent with the 
relation {R a bcd R abcd ) 1 / 2 > 4\^2m/r 3 derived by Christodoulou for the case of the collapse 
of a massless scalar field in [9]. 

Interestingly, in spite of the divergence of the Ricci scalar, the quantity aj3 r 2 R sr remains 
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Figure 9: Constraint preservation on a time level surface near the singularity. Left: Misner-Sharp mass 
(narrow line), the absolute value of its p-derivative (dotted line) and the numerical violation of constraint 
(|26[) (thick solid line). Right: the a metric component and the numerical violation of constraint (1791) . 

finite. Thereby the Einstein-Hilbert action 

S EH = R sc e = ^\ / a[3r 2 R sc dpdT, (88) 

remains finite as well, for any particular choice of pi,P2 with < p\ < p 2 < p ma x m 
S' = [pi,p 2 ] xS 2 . 

7.4 Gravitational collapse with spatial topology S = § 3 

Whenever the topology of the initial data surface, So, is chosen to be S 3 = [0, tt] x S 2 , there 
has to be two origins on So in the considered spherically symmetric setting. Motivated by 
and mimicking the simplicity of the geometric setup of the FRW cosmological model, we 
used procedure B (see section I7.ip and chose the functional form of the initial a and r as 

a = Rq and r = Rq sin p, (89) 

To have a nearly central symmetric time evolution, the initial data for the scalar field was 
chosen to be slightly asymmetric by requiring ip and ip T to take the form 

ip = Vo cos np + tpi cos n'p and ijj T = ipo cos np + ipi cos n'p , (90) 

The parameters i?o, tpo, V'l) V'Oj ipi, n and n' in the above functional expressions were fixed 
as 

#0 = 10, = 0.4, tpx = 0.005, ip = 1.2, ipx = 0.01, n = 3 and ri = 2, (91) 

while the number of the spatial grid points N was chosen to be 1,000. 

On Fig.[2]the time evolution of the gravity-matter energy density distribution associated 
with a shell of radius p, £ GM , is shown. The left panel of Fig. [15] shows a somewhat exotic 
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Figure 10: Constraint preservation in case of dynamical lapse. 



distribution of the marginally trapped surfaces produced by the time evolution. This plot 
justifies that the formation of the curvature singularity is censored, i.e., a connected future 
trapped region develops in advance to the appearance of the singularity. There are two 
apparently smooth curves starting at the bottom and ending at the opposite upper corner 
where and 0+ are identically zero, respectively. These curves intersect at p = tt/2. The 
corresponding 2-sphere is a maximal surface. The points on the thick part of the curves 
represent future marginally trapped 2-surfaces while the thin parts depict the location of 
the past marginally trapped 2-surfaces. Thereby the associated past apparent horizon may 
also be referred as the boundary of the dynamical white hole region. 

On the right panel of Fig.rTS]the time dependence of the spatial 3- volume, V = Jq y/afi&p, 
is shown. It increases exponentially up to r ~ 0.48 after which the tendency is reversed and 
an even faster decay of V is experienced. 

An interesting new feature of the evolution is that the metric function a appears to 
blow up while approaching the singularity. However, its blow-up is compensated by the 
decay of /3; the proper time t(p) = /J* ^fafi&T was found to be finite along any constant 
p world-line — meaning that the singularity is in finite proper time distance from the initial 
data surface. 

The r-dependence of the square root of the Kretschmann scalar and the Misner-Sharp 
mass, along the world-line p = 2.5, are indicated on the left panel of Fig.[16j On the right 
the p-dependence of the critical exponents 7x, 7m and 7^ — 7 m are given. Interestingly, 
Christodoulou's relation {R a bcd,R abcd ) 1 ^ 2 > 4y/2m/r 3 — derived only for the case of the 
collapse of a massless scalar field in a spacetime with spatial topolog y R 3 in [9] — does apply 
here as well. 

The value of the Einstein-Hilbert action, 



remains finite as it is indicated by the r-dependence of the Lagrangian L = Air a(3 r 2 R sc dp, 




(92) 



25 



Figure 11: Time evolution of the gravity-matter energy density distribution, S GU . The black curve 
indicates the apparent horizon. [Coloured and interactive versions of this and some other 3-D figures can be 
find at http://www.kfki.hu/~cspeter/numrel/2009-ekg/index.html ] 

shown on the left panel of Fig. 1151 

7.5 Gravitational collapse with spatial topology S = S 1 x § 2 

Whenever the topology of the initial data surface So is S 1 x § 2 , the simplest initial configu- 
ration to start with is a homogeneous 'torus' with constant initial radius and with constant 
initial expansion. Accordingly, we used scheme A and chose the initial r and r T as 

r(p) = r and r T (p) = r , (93) 

while the number of the spatial grid points N was chosen to be 4,000. 

For the the initial ip and ip T , we used the relations in (|90p with the parameter values in 
(|91[) . while the values of the metric parameters were chosen as 

r = 1, and r = 1 . (94) 

It is important to keep in mind that now p is a periodic coordinate, along the S 1 factor. 
This periodicity length was chosen to be 2ir. Notice also that according to the above choice 
of initial data, there is no origin at all on So- 

On Fig.[T7]the time evolution of the gravity-matter energy density distribution associated 
with a shell of radius p, d? GM , is shown. As opposed to the other cases, S GM tends to 
zero at the singularity. The Kretschmann scalar blows up in the singularity such that 
Christodoulou's relation {R a i )Cl iR abcd ) 1 ^ 2 > 4\^2m/r 3 holds. The location of the future and 
past apparent horizons are also indicated on the left panel of Fig.Q2] . 

The appearance of curvature singularities are preceeded by the formation of trapped 
regions. As in the previous case, both future and past apparent horizons are formed. They 
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Figure 12: On the left, the r — singularity (dashed line) and the apparent horizon (continuous line) are 
shown. On the right, the p-dependence of the Misner-Sharp mass along this apparent horizon is plotted. 

intersect at the 2-surfaces with p ~ 0, n/2, tt, 3tt/2. Interestingly enough, while the 2-surface 
with p ~ and p ~ ir are maximal the other two 2-surfaces have no definite character. 
More precisely, they are locally maximal in the spatial p-direction whereas they are locally 
minimal in the timelike r-direction. The spacetime region with untrapped surfaces is very 
limited in the present case. It is represented by the points between the future (thick solid 
line) and the past (thin solid line) apparent horizons on the left panel of Fig. 1181 

The r-dependence of the Lagrangian L = 4ir aj3 r 2 R sc dp , indicates (see the right 
panel of Fig. ll8[) that the Einstein-Hilbert action, 



remains finite as in the previous cases. 

Finally, the r-dependence of the spatial 3-volume, V = Jq V yfafldp, is also shown on 
the right. The volume starts to increase rapidly, then a short oscillation followed by an 
extremely fast contraction can be seen. 

8 Summary 

In most of the former numerical simulations, one of the aims was to minimise the extent 
of the trapped region. This was achieved by making use of techniques like singularity (and 
trapped region) avoiding slicings and black hole excision — the latter was originally suggested 
by Unruh in 1984 (for its first adaptation in numerical simulations see, e.g., [431 139] ). The 
unsatisfactory aspect of this trouble avoiding attitude, and the reason for choosing here the 
simplest possible framework of spherically symmetric dynamical configurations are justified 
by the following comments of David Hilbert (1902): 

"In dealing with mathematical problems, specialisation plays, as I believe, a still more 
important part than generalisation. Perhaps in most cases where we unsuccessfully seek the 




(95) 
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Figure 13: Left: Blow up rates of the square root of the Kretschmann scalar, R a bcd R abcd , the Ricci scalar 
curvature R sc = g ab R a b, and the Misner-Sharp mass m along a timelike curve, the p — 2.5 line. Right: 
p-dependence of the critical exponents of the blow-up. 

answer to a question, the cause of the failure lies in the fact that problems simpler and easier 
than the one in hand have been either incompletely solved, or not solved at all. Everything 
depends then, on finding those easier problems and on solving them by means of devices as 
perfect as possible and of concepts capable of generalisations. " 

The key technical achievements that, as we believe, have not been applied before in 
numerical simulations, are the following: 

(1) A strongly hyperbolic (symmetrisable) first order system of evolution equations was 
singled out for 4-dimensional spherically symmetric gravitating systems. 

(2) The analytic setup ensures that time evolution can be studied on equal footing in 
trapped and untrapped regions. 

(3) The numerical framework applies this analytic setup and incorporates the techniques 



By making use of these technical developments we achieved the following results: 

(1) By introducing a suitable evolution equation for the lapse function (3, the extent of the 
investigated spacetime domain was enlarged significantly and the physical singularities 
were approached arbitrarily closely everywhere. 

(2) The location of the future and past apparent horizons has been determined. 

(3) Detailed investigation of the rate of curvature blow-up while approaching the singu- 
larity. 

(4) The Einstein-Hilbert action remained finite in all the investigated cases, in spite of 
the blow up of the Ricci scalar. 



of AMR. 
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Figure 14: The time evolution of £ GM . The black and white curves indicate the location of the future and 
past apparent horizons. 

As one of our aims was to clear up some conceptual issues of topology change, let us 
close this paper by related comments. 

In dealing with the problem of topology changes, one should start by recalling the result 
of Geroch [2jj which asserts that if topology change develops then there must exist either 
closed causal curves beyond the Cauchy horizon^ or a spacetime singularity has to appear. 
Due to Tipler's theorem [52], the latter case manifests itself if Einstein's equations are 
imposed. Based on these results, it has been widely held that no indication of topology 
changes will ever show up in classical general relativity, therefore the quantum theory should 
be investigated to see whether they may occur. 

The idea that the topology may change in quantum gravity was originally proposed by 
Wheeler [451 146] . Since then, one of the most important questions in any quantum theory 
of gravity is whether there is a non-zero probability for the topology of space to change. 
There have been a number of classical investigations aiming to demonstrate the feasibility of 
topology change by making use of the techniques of differential geometry and topology (see, 
e.g., jU EH] [T71 [4jJ ) . Nevertheless, no quantitative investigations have been carried out yet. 
Therefore it is important to emphasise that some of our findings provide the first definite 
quantitative support to all the former speculations concerning the existence of topology 
changes. 

In summarising our pertinent observations, we can say, in accordance with the above 
recalled results of Geroch and Tipler, that instead of having a regular Cauchy horizon, 
along with a causality violating region beyond, a spacetime singularity develops at the 

9 As it was discussed in the introduction in the present context one could always think of the limit of the 
time level surfaces as the Cauchy horizon. 
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Figure 15: Left: The location of the future (thick solid line) and the past (thin solid line) apparent horizons, 
along with the singularity (dashed line). Right: The time dependence of the 3- volume, V = yfaP&p, and 
the Lagrangian, L — 4n a/3r 2 R sc dp. 

"new origins". It was also found that the Kretschmann scalar always blows up there. It is 
important to emphasise that the Einstein-Hubert action S EH = Jm=Ex[to T j R sc e was found 
to be bounded in all of our investigations. 

The existence of time developments with apparent topology changes may be significant 
in quantum theoretical considerations. For instance, in the sum over histories approach 
to quantum gravity, the transition amplitude for the topology change between the two 
Riemannian manifolds (Ei,^i) and (S2,/i2) is supposed to be given by the formula 

((i: 1 ,h 1 )\(E 2 ,h2)) =J2 [ exp[iS EH ]Vg ab , (96) 

where the boundary of M is the disjoint union of Si and S2, S EH denote the Einstein- 
Hilbert action and is the space of 4-dimensional Lorentzian geometries on M, while 
the sum is over all 4-manifolds whose boundary is the disjoint union of Si and S2 with 
Riemannian metrics h\ and /12 (see, e.g., [4"5 | |4*6 ] \27\ l28j [18]). 

Admittedly, the right hand side of (|96p is completely formal and is far from being defined 
as yet. Nevertheless, regardless of the specific form of the measure on the space of metrics, 
T>g a t,, it is widely held that not only the smooth Lorentzian metrics but all the metrics 
with finite Einstein-Hilbert action should be taken into account in evaluating the functional 
integral. Therefore the spacetimes with topology changes investigated in this paper, might 
be of interest in quantum theoretical considerations. 

It is of obvious importance to know whether the methods introduced in this paper could 
be adopted in more generic geometrical setup allowing the presence of gravitational waves 
and with the inclusion of binary black holes. We would like to mention that by making 
use of the analytic setup suggested in [31] such a generalisation seems to be possible. The 
results of our corresponding investigations will be published elsewhere. 
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Figure 16: Left: The r-dependence of (Rabcd R abcd ) 1/2 and the Misner-Sharp mass m along a timelike 
curve, the p — 2.5 line. Right: p-dependence of the critical exponents. 
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